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ABSTRACT 


The finite element method is used to investigate the static behav- 
ior of laminated composite flat plates and cylindrical sheila. The 
structures are composed of an arbitrary number of layers with arbitrary 
material properties and fiber direction within each layer. The analy- 
sis incorporates the effects of traverse shear deformation in each 
layer through the assumption that Che normals to the undeformed layer 
midsurface remain straight but need not be normal to the mid-surface 
after deformation. 

The present finite element formulation is based on assumed first 
order polynomial displacement functions. Itiie nodal displacements are 
obtained by minimizing the corresponding potential energy to establish 
approximate equilibrium equations. The solution of these simultaneous 
equations is then the ultimate goal for the study. A digital computer 
program which utilizes part of the available SAP V computer program was 
developed to perform the required computations. The program includes 

s 

a very efficient equation solution code which permits the analysis of 
large size problems . 

To assess the behavior of the developed finite element^ solutions 
were obtained for some laminated composite plate and cylindrical shell 
problems for which analytic solutions are available. The accuracy of 
the stress and displacement results obtained by using the developed 
finite element is determined for different plate and shell proportions . 
The rate of convergence and limitations on the use of the developed 
element are also discussed. The method is finally applied to the prob- 
lem of stretching and bending of a perforated curved plate. 
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CHAPTER I 


INTRODUCTION 


Thin shells have become a popular subject In the field of struc-* 
tural mechanics In the past two decades. Numerous books and papers on 
their analysis have been published because of their use In such struc- 
tures as aerospace structures, pressure and underwater vessels, and 
nuclear reactor structures subjected to both static and transient laods. 
The advent of advanced flber-relnforced composite materials which has 
been called "the biggest technical revolution since the jet engine" (1), 
has led to increasing use of lamlnted shells. 

In flber-relnforced composites such as boron-epoxy and graphite- 
epoxy combinations, the fibers provide the majority of the strength and 
stiffness. The function of the matrix Is to support and protect the 
fibers as well as to distribute and transmit load between the fibers. 

The latter function Is especially Important If a fiber breaks, for then 
load from one side of the broken fiber Is transferred to the matrix and 
subsequently to the other side of the broken fiber and to the adjacent 
fibers. The matrix Is of such low shearing stiffness compared to the 
fibers, however, that shearing deformations become important in appli- 
cations involving loads normal to the structure. 

In laminated structures the desirability of accounting for shear 
deformations in each layer has made it necessary to extend the classical 
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analyses based upon the Klrchho££-Love hypothesis which neglects the 
e££ects o£ transverse shearing de£ormatlon entirely (re£erences 3,4) or 
upon the Tlmoshenko-Mlndlln modl£lcatlon which deals with an average 
shearing deformation (references 5-11) . The Importance o£ a more rigor- 
ous analysis is indicated by the work of references 12-15 which investi- 
gates particular elasticity problems of the layered plates by means of 
exact analytic solutions. 

For more general problems It Is difficult to obtain an analytic 
solution, however, exact or approximate. Due to developments in compu- 
ter science and digital computers, the finite element method offers the 
possibility of yielding numerical results for problems unsolvable by 
other methods. Some investigations utilizing the fi n ite element method 
for laminated plates including shear deformations are already available 
in the literature. In references 16 and 17 the finite eleme n t model 

t 

includes only average shearing deformation effects. Mau, Tong, and 
Plan (reference 18) provide for general laminate behavior by allowing 
each layer to undergo independent rotations of the normal to its unde- 
formed middle surface in a hybrid finite element model. A modified hy- 
brid finite element has been developed by Spilker, Chow, and Orringer 
(reference 19) who seek to obtain reasonable accuracy with a consider- 
ably reduced number of nodal degrees of freedom. 

Because the element used is generally a rectangular plate element 
in these investigations, however, the solution of such problems as lami- 
nated plates with perforations for the determination of stress concen- 
tration effects around the hole cannot be readily studied. In the 
present investigation the stiffness matrix for a triangular curved 
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element which has better application to these and other problems is 
derived and applied to various problems. The study Is limited to a 
cylindrical shell and flat plate element which consists of a right 
triangle In the lines of curvature coordinate system and with the per- 
pendicular sides oriented in the directions of the lines of curvature 
coordinate curves. For other shells and element shapes a numerical 
Integration technique (reference 20) Is possible. 

In the study a displacement finite element model rather than the 
similar hybrid stress model of previous Investigations is used . The 
motivation for this simplified model was a belief that much the same 
accuracy as the hybrid stress model could be obtained, with a consider- 
able simplification of the calculations involved. This belief is borne 
out to some extent by a comparison of displacements and stresses for 
several problems solved numerically by the two finite element models, 
in that much the same results can be obtained, at the price, however, of 
a slower rate of convergence. 

The finite element model is based upon a theory of layered plates 
and shells which retains the Love-Klrchhoff assumption of a normal 
deformation which is not a function of position in the shell thickness 
coordinate. The assumption that straight lines normal to the undeformed 
shell remain normal to the deformed shell is replaced by the assumption 
that straight lines in each layer normal to the unformed shell remain 
straight but rotate with respect to the normal to the deformed shell by 
an amount which varies from layer to layer. A similar assumption has 
been utilized for beams by Swift and Heller (references 21-24) and by 
Durocher and Solecki (references 25-26) for two and three layer plates. 


The limiting case of an Infinite number of Infinitesimally thin layers 
has been shown to yield good agreement with the results of three 
dimensional isotropic plate theory by Boal and Relssner (reference 27) . 
For a single layer plate the theory reduces to a form similar to that 
of Llbove and Batdorf (reference 28) for sandwich plates. 

The derivations that lead to the finite element stiffness formula- 
tion are discussed in Chapter II. A computer program based upon the 
results for circular cylindrical shells and flat plates is discussed in 
Chapter III and listed in Appendices A and B. Finally- the results for 
several problems investigated with the present analysis are discussed 
in Chapter IV. 
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II. THEORETICAL BACKGROUND 

2«1 Assumptions 

The thin shell theory adopted here is associated with the following 
set of assumptions: 

1. The shell is composed of layers of orthotropic material, the 
principal material directions of which do not necessarily coincide with 
the coordinate directions. (Fig. 2-1). 

2. There is a linear distribution of tangential displacements 
t&ough the thickness of each layer. The linear distribution may vary 
from layer to layer. (Fig. 2-2). 

3. The normal displacement is constant through the entire thick- 
ness of the shell at each point on the shell reference surface, provided 
there are no portions of one or more interior layers missing. 

The various elements of the theory are as follows. 


2.2 Strain-Displacement Relations 

The strain-displacement relations used are obtained from three- 
dimensional elasticity theory. A system of orthogonal curvilinear coor- 
dinates for points within the shell wall is defined by the coordinates 
a and 3 corresponding to the lines of curvature on the bottom surface 
and distance ^ along the normal to the surface . (Fig. 2-3). Then, from 
reference 29, the following equations can be established: 
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where u and Uo are tangential displacements, w is the normal displace- 

ct P 

ment, g and g^ are first fiindamental quantities for the lines of curva- 
a p 

ture coordinate curves on the reference surface and R^ and R„ are normal 

cc p 

cuirvatures of the reference surface in the directions of the coordinate 
curves. 

The basic assumptions in this study are those of linear tangential 
displacement and constant normal displacement distributions through the 
layer thickness. For the ith layer, then. 
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where u and Uq are tangential displacements of the kth Interface and 
a p 

are functions of a and 3 only. The coordinate C Is now measured In 
each layer from the bottom of the layer so that 0 < C < h^. Eqs. (2.1) 
then become 
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2.3 Material Properties 

The material of each layer is assumed to be anisotropic in surfaces 
parallel to the reference surface of the shell. In accordance with the 
usual assumptions of plate and shell theory, layer strains are assumed 
independent of transverse stresses. 

The stress-strain relations are then taken as 

{a} « [E]{e} C2.5a) 

where {or} and {e} are respectively the stress and strain matrices 
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and E Is the symnetrlc material constant matrix 
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Poisson's ratios, and shear moduli as discussed in detail in Ref. (2). 
For an orthotropic material with the principal material directions coin- 
ciding with the coordinate axes. Equations (2.5a) may be written as 
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The principal material directions for an orthotropic material often do 
not coincide with the coordinate directions, however, especially in a 
laminated composite structure whete the orientation of material can be 
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varied from layer to layer. Therefore it is convenient to effect a 
transformation of the stress-strain relations from local material axes 
to the overall coordinate axes. Let the angle between the global and 
local orthogonal layer material axes be denoted as 0 (Fig. 2-4). Then 
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These transformations are derived in most texts on strength of materials 
and theory of elastiicty, for example Ref. 

Then 
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2.4 Principle of Minimum Potential Energy 

The theorem of virtual work states that if an elastic body is in 
equilibrium under the action of prescribed body and surface forces, the 
work done by these forces in a small additional displacement is equal to 
the change in internal strain energy. From the theory of elasticity the 
internal strain energy for a layered shell is given by 
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With the stress-strain relations given by Eqs. (2.5) the strain 
energy may be written as 
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where the surface integral is over the area of the reference surface and 
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The components of the strain matrix for each layer are given by Eqs. 
(2.3) and [E^] is the material modulus matrix for each layer of the form 
of Eq. (2.5d). Multiplication of the matrices in the integral then 
yields the following expression for the strain energy 
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If the applied forces (Fig. 2-5) are conservative the work done by 
these forces is given by 
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where 1/R is the normal curvature of the reference in the direction of 
3 

the boundary curve and is the twist of the surface referred to 

surface coordinates normal and tangential to the boundary curve, while 
^ns* ^n? applied stresses parallel, normal, and tranverse to 

the boundary cuxrve, respectively. 

The principle of virtual work can now be stated as follows as the 
principle of minimum potential energy: 

Of all displacements satisfying the given boundary condition on dis- 
placement, those which satisfy the equilibrium equations make the poten- 
tial energy 11^ ■ U-W stationary. An expression such as 11^ is called a 
functional to imply that is not simply a function of displacements 
and displacement derivatives but depends on their integrated effect. 


2.5 The Finite Element Method 

The principle of minimum potential energy is the basis for some 
procedures for obtaining approximate solutions for problems which cannot 
be solved analytically. The older classical procedure kno^ra as the 
Rayleigh-Ritz method (ref. 31) expresses the solution in terms of.series of 
known trial functions which each satisfy geometric boundary conditions 
^hich are continuous and have continuous derivatives over the entire 
region of interest. The coefficients which multiply these functions in 
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the series are determined by the condition that the potec.tial energy be 
a minimum, yielding the conditions that the derivatives of the functional 
Up with respect to each of the coefficients must vanish. 

The finite-element method is a variation of this procedure which 
seeks to define the solution in terms of different trial fxinctions, 
usually polynomial expressions, for different discrete regions of the 
structure. Only the trial fxinctions and, in some cases, certain of 
their derivatives are continuous at the common boundary of adjacent 
regions. Such a series of functions is said to be a piecewise continuous 
fit to the solution. The solution within each discrete region is ex- 
pressed in terms of displacements and, if necessary, derivatives at 
discrete points or nodes, usually on the boundary. Thus while the form 
of the polynomial expression is the same for each region, the coeffi- 
cients vary from region to region. The potential energy of the entire 
structure is then obtained in terms of these nodal quantities which take 
the place of the function coefficients in the classical Rayleigh-Ritz 
procedure. Extremallzation of the potential energy then yields equations 
for determination of the nodal quantities . (for example refs. 32 and 33). 

Tlie advantage of the finite element method over the classical 
Rayleigh-Ritz method is that the same discrete region functions are 
valid for any boundary conditions on the structure. Thus the algebraic 
manipulations yielding the form of the potential energy expression need 
be carried out only once, after which an algorithm to yield the results 
for regions having the same general shape but different parameters may 
be written for automated digital computer applications. 


A 


t 





2.6 Finite Element Formulation 

The finite element model usually chosen to represent a plate or 
shell structure is a series of triangles or rectangles since the analysis 
of the structure is reduced to that of a two-dimensional continuum. For 
the present problem the inclusion of shearing effects in each layer re- 
quires the analysis of a number of interacting two-dimensional contlnua, 
one for each layer interface. Thus an element consists not of a single 
rectangle or triangle but a number of these, one above the other, for 
corresponding regions on the layer interfaces. While a rectangular ele- 
ment shape is adeqtsate for a structure with regular boundaries, the 
shape of the element is chosen to be a triangle in the a, 3 plane to be 
better adaptable for the analysis of regions with irregular boundaries, 
as shown in Fig. 2-1. 

If the vertices of the triangle element are numbered as shown in 
Fig. 2 - 73 . and the shape of the triangle is such that and 

the corresponding strain energy of the triangular element of the 
1^^ layer is given by 



i J j 

If j however, the vertices of the triangle are numbered so that 

but 33<S2 (Fig. 2-7b), the limits of integration on 3 are reversed. 

In evaluating the strain energy of a typical triangular layer ele- 

til 

ment the displacements in the j element of interface i are assumed 
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to vary linearly within each element as 


“a^ - 

+ b«6 + 

wM - 33^0 + b 3 ^S + C3^ 

The coefficients may be obtained in terms of the comer 
ment by means of the equations 
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(2.16) 


displace- 


(2.17) 


(2.18a) 


(2.18b) 


(2.18c) 
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In addition the numbering system shovm. in Fig. 2-3 is introduced. Then 


a relationship between the element strain matrix {i^^jand the displace- 
ment matrix 

M' 


idl 


'ij 


u 


ij 
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(2.19a) 


can be obtained as 




(2.19b) 


where 
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and the quantities Rj^ are the components of the matrix 


(B2-S3)a + (03-02)6 + OjBj-o^Bj (63-6j)o + (o}-o^)6 + o^bJ-oJBj (6j-62>o + (a^-oJ)B + oJbJ- o^Bj 
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(2.19c) 
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Note that only the last two rows of the matrix [G] change from layer to 
layer. Then the layer element strain energy is given by 


{dLCKl.. {d> 


'ij ^“^ij 


(2.20a) 


with 
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S)]dade 

ij 


and 


[F(ot,3)]y - 


(2.20b) 

(2.20c) 


The element stiffness matrix may be obtained in an integrated 


form if the components of the first and second fundamental forms have a 
simple variation with the lines of curvature coordinates a and 3 and if 
one has the patience to carry out the matrix multiplications and the 
integrations. An alternative approach is a computer program which will 
carry out the matrix multiplications at selected points and perform the 
required integrations numerically. 

Similar operations will yield the work done by applied surface 
loads. For example, work of normal surface loading on any element is 



ftj— ftj 
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(2.21a) 



where 
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(2.21b) 


and Wj is given by 


w^ - W2J { 
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(2.21c) 


If is assumed to vary linearly within the element as an approximation 


to the actual distribution, is given by an expression similar to that 


of Eq. (2.21c). Then 
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The potential energy is now given by 
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E E - E AW, 
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(2.22a) 


(2.22b) 


V 


(2.22c 


and Is expressed In terms of the displacements at the nodal points of 
the structure. Minimization of the potential energy with respect to 
each of these nodal point displacements then yields a set of linear 
equations for their determination. Convergence of the process to the 
correct solution occurs as the size of the elements become smaller with 
an Increasing number of elements. 

With the displacements known at the vertices of each triangular 
prism layer element, the stress distribution within each element may 
be obtained from the relationship 
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A better approximation to the transverse shear stresses x and x„ can 

ap 3 p 

be obtained by using the equations of equilibrium of stresses within 
each layer as In reference 29. 


2.7 Specific Results for a Right Triangle Element of a Circular 
Cylindrical Shell and a Flat Plate 

For a right triangle element of a cylindrical shell, with a taken 
as the coordinate x along the generator and 0as the clrciimferentlal 
angle 9 so that 


the matrix is given by 


1 

R 

OO 

R 
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In addition 
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IHie stiffness matrix is given by 
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while the work expression of Eq. (2.22b) is given by 



r ^ 



1/12 

1/24 

1/24 



‘^2j 

1/24 

1/12 

1/24 

< 



1/24 

1/24 

1/12 




RA0^Ax^ 


(2.28) 


The Integrations of Eq. (2.27) have been carried out and have been 
Incorportated In the computer program listed In Appendix A. In doing 
so care was taken to distinguish between the four differently oriented 
right triangles shown In Fig. 2-9 for which various terms in the com- 
ponents of the stiffness matrix change sign. 

Similar results hold for a flat plate for which R is infinite 
while R A0 » Ay . Then 


(2.29a) 

(2.29b) 


l2j^(x,0) - l2^(0,x) » -|h^ 
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l2^(x,0) - l2^(0,x) - l^^(x,0) » l^^(0,x) * -|h^ 


(2.29c) 
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with 
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(2.31b) 
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(2.31c) 

«12 

■ Ay^ (x-x|)-Ax^ (y-y|) 

(2.31d) 
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(2.31e) 


and the stiffness matrix is given by 
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^ Ax^ 


[K] 


ij 
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[G] [B] [G] dxdy 


(2.32) 


while the work expression is identical with Eq. (2.28) with RA0^ 
replaced by Ay^ . 


CHAPTER III 


COMPUTER PROGRAM FOR THE ANALYSIS OF LAYERED PLATES AND SHELLS 

The finite element theory of the preceeding chapter using the 
particular right triangular element developed for flat plates and 
cylindrical shells has been incorporated into the computer code listed 
in Appendix A for which input and output data are given.in Appendix B. The 
code utilizes a ntimber of subroutines taken from the SAP I V and SAP V 
codes (references 23-25) together with some additional subroutines de- 
vised for the present program. The SAP system codes are particularly 
efficient subroutines for the processing of input data, the generation 
of the system of equations to be solved, and the solution of an equa- 
tion system of large size and bandwidth. The added subroutines deal 
with the input and output peculiar to laminated plates and shells. 

The subroutines and their functions are as follows: 

1. HEDPRT - Reads the total number of nodal points, 

elements, layers, and the type of analysis. 

2. INPUTJ - Reads or generates the coordinates and the 

restraint conditions of each nodal point, reads 
the material properties and thickness of each 
layer, generates the variable number of each 
degree of freedom. 

3. LAYSHL, AYSHL - Reads or generates the layer element type 

and the nodal points forming each element, cal- 
culates the element stiffness matrix, calls ADDSTF 
and CALBAN 
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4. 


INL - Reads or generates nodal forces. 

5. CHECK - Verifies sufficiency of input data. Prints an 

error message if criteria not met. 

6. ADDSTF- Stores layer element stiffness matrix in proper 

location in overall stiffness matrix, divides 

the stiffness matrix and load vector into blocks 

for storage. 

7. CALBAN - Calculates bandwidth of system of equations. 

8. SOLEQ - Controls input and output data from subroutine 

SESOL. Calls PRIND, STRESS, STRSC. 

9. SESOL - Solves positive definite symmetrical system 

equations by Gaussian elimination. 

10. ERROR - Prints error message if limit on storage require- 

ments is exceeded. 

11. STRESS, STRSC - Calculates stresses for each element. 

12. PRIND - Prints out all required input and output. 

The input data required for the static analysis of layered plate 
and shell structures is given t~i detail below. The nodal points should 
be numbered consecutively along straight, normal lines from the bottom 
layer to the top layer before proceeding to another nodal point on the 
bottom layer. The elements should! also be numbered consecutively 
through the thickness before proceeding to another element on the bot- 
tom layer. The computations were carried out on the USC School of 
Engineering DEC- 10 System Computer. 
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Heading Card 

(12A6) 


Columns 

Variable 

Entry- 

1-72 

HED(12) 

Enter the heading information to be 
printed with the output 

Master Control Card (415 

, 14, 11, 415, 2X, 311, 315, 5X, PlO.O) 

Columns 

Variable 

Entry 

1-5 

NXJMNP 

Total number of nodal points in the model 

6-10 

NELTYP 

1 

11-15 

LL 

Number of Structure load cases 

16-20 

N F 

0 

21-24 

KDYN 

0 

25 

NDYN 

0 

26-30 

MODEX 

Program execution mode 

0 if problem solution wanted 

1 if only data check wanted 

31-35 

NAD 

0 

36-40 

KEQB 

0 

41-45 

NIOSV 

0 

48 


0 

49 

GENFRT 

0 if printout of overall stiffness 

matrix not wanted 

1 if printout of overall stiffness 

matrix wanted 
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Columns 


1&rlable 


Eutiry 



50 ELPRT 0 if printout of element stiffness 

matrices not wanted 
1 if printout of element stiffness 
matrices wanted 

51-55 IWTC 0 

56-60 MINBND 0 if bandwidth, minimization not wanted 

1 if bandwidth minimization wanted 
61-65 IPLT 0 

66-70 NRSC 0 

71-80 GRAV 0 

III. .Structural Data (415) 

Columns Variable Entry 

1~5 NELTYP 13 if cylindrical shell element 

14 if flat plate element 
6-10 NOME Total number of elements 

11-15 NUML Total number of layers 

16-20 NUMPX Total number of restrained nodal points 

IV. Material iroperties (P5.1, F3.1, 9F8.0) (One for each Layer) 

I 

Columns Variable Entry 

1-5 A Angle of the layer material axes of 

orthotropy to global cnordlnate 
system (Equal to 0.0 if material 
is generally anisotropic) 
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Columns 


Variable 


Entry 


6-8 

H 

9-16 

E 

11 

17-24 

^2 

25-32 ■ 

®13 

33-40 

^22 

41-45 

®23 

49-56 

E 

33 

57-64 

) 

G 

44 

65-72 

G 

45 

73-80 

G 

55 


Thickness of the layer 


Material elastic moduli (E, . ■ E 

13 23 

» G ■ 0.0 if material is ortho- 
45 

tropic and properties are referred 
to axes of ortho tropy) . 


Nodal Point Data (13, 3F7.4) (One for each nodal point) 

Coltimn Thriable Entry 

1-3 I Assigned number of the nodal point 

numbered consecutively through 
the structure thickness before 
proceeding to next node in bot- 
tom face 


4-10 

X(I) 

X -ordinate 

11-17 

Y(I) 

Y-ordiante 

1&-24 

Z(I) 

Z-ordinate 
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\lAj-J . . i . . . : . J : J ‘ ' ' ! . 

d 

G 

VI. Restrained Point Data (13, 311) (One for each restrained nodal point) 


U Columns 

Variable 

Entry 


N 

Number of restrained nodal points 

li 4 

ID(N.l) 

Restraint condition in x- translation 

n 


1 if point is restrained 

rn 


0 if point is not .restrained 

Q 5 

ID(N,2) 

Restraint condition in y-translation 

n 


1 if point is restrained 

u 


0 if point is not restrained 

0 

ID(N,3) 

Restraint condition in a-translation 

n 


1 if point is restrained 

0 


0 if point is not restrained 

J VII. Radius For 

Shell Element 

(F13.3) (outfitted for flat plate 


analysis) 

VIII. Element Load Mxiltlpliers (4F10.0) 


Columns 

Variable 

Entry 

1-10 

EM(1) 

0.0 

11-20 

EM(2) 

0.0 

21-30 

EM(3) 

0.0 

31-40 

EM(4) 

0.0 


LJ 


U 

n 

l1 
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Element 

Data (415,311) 

(One for each element) 


Columns 

Variable 

Entry 

Notes 

1-5 

M 

Element Number 

(1) 

6-10 

11 

Vertex of the element 

(2) 

11-15 

Jl 

Vertex of the element 


16-20 

K1 

Vertex of the element 


21 

L 

Layer number of the element 

(2) 

22 

MCD 

Element Code 

(2) 

23 

ML 

0 without load on the element 

1 with load on the element 

(3) 

(4) 


Notes: (1) Numbered consecutively through the structure thickness 

before proceeding to next bottom layer element 

(2) There are 4 different elements based on their orientation 
as shown in Fig. 2-9. In all cases 


12 - II + 1 
J2 - J1 + 1 
K2 - K1 + 1 

(3) When there is no load acting on the element, input 
the next element data 

(4) When there is any load acting on the element, input the 
load data (next input) before the following element data 


X. Load Data (18F4.2) (One for each loaded element) 

Columns Variable Entry 

4- 7 QX(Il) Applied x-direction stress at II 

5- 8 QY(Il) Applied y-direction stress at II 

Applied x-direction stress at Jl 


9-12 


QX(Il) 

QY(Il) 

QX(Jl) 


.,11 ^ — li ■ II ii tfrfL i 




Coliimns 

Variable 

Entry 



13-16 

QY(J2) 

Applied y-direction stress 

at J2 


17-20 

QX(Kl) 

Applied x-dlrection stress 

at K1 


21-24 

QY(Kl) 

Applied y-direction stress 

at K1 


25-28 

QX(I2) 

Applied x-dlrection stress 

at 12 


29-32 

QY(I2) 

Applied y-direction stress 

at 12 


33-36 

QX(J2) 

Applied x-dlrection stress 

at J2 


37-40 

QY(J2) 

Applied y-direction stress 

at J2 


41-44 

QX(K2) 

Applied x-direction stress 

at K2 


45-48 

QY(K2) 

Applied y-direction stress 

at K2 


49-52 

QZ(Il) 

Sun of applied z-direction 

stresses 

at 



nodes 11 and 12 



53-56 

FZ(Ia) 

Sum of applied z-directions 

concentrated 



forces at nodes 11 and I2 



57-60 

QZ(Jl) 

Sum of applied z-direction 

stresses 

at 



nodes J1 and J2 



61-64 

FZ(Jl) 

Sum of applied z-directions 

concentrated 



forces at nodes J1 and J2 



65-68 

QZ(Kl) 

Sum of applied z-direction 

stresses 

at 



nodes K1 and K2 



69-72 

FZ(Kl) 

Sum of applied z-direction 

concentrated 



forces and K1 and K2 
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XI. Load Point (13) (One for each loaded point in every load case 


could be in any direction) 


t I TiiiritWw 






IV. FINITE ELEMENT ANALYSIS OF SPECIFIC PROBLEMS 

A number of problems have been InvestigaCed, some Co assess the 
accuracy and convergence properties of the finite element analysis 
derived and others to explore the usefulness of the analysis for prob- 
lems difficult to solve by other methods. A description of the various 
problems and the results obtained are as follows: 

4.1 Cylindrical Bending of a Wide Simply Supported Laminated Plate 

Subjected to Sinusoidal Surface Pressure 

Pagano in Ref. 13 has given an exact plane strain elaiJticity 
solution for Che problem of cylindrical Bending of a layered plate sub- 
jected to normal pressure. (Fig. 4-1) . The material of each layer is ortho tro- 
pic and identical. The directions of the axes of orthotropy of alternating 
layers vary by 90®, however. The boundary conditions are those of 
vanishing edge normal displacement and stress. For a surface loading 
which varies sinusoidally, then, the displacements and stresses have 
sinusoidal variations as well. 

To obtain an exact solution of the problem within the framework of 
Che present theory the transverse displacement is taken in the form 

w ■ W sin (4.1a) 

while the longitudinal displacements are taken in the form 

u^ » cos ^ (4.1b) 

XX li 
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Then Eq. (2.12c) yields 


..i IT , ITX 

-U — sin ~ 
X L L 



0 

0 

0 





cos 


TTX 

L 


0 


(4.2) 


The use of Eqs. (2.12), (2.14) and Eq. (3.2) results in the follow- 
ing expression for the potential energy per unit width of the layered 
plate 
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where is the amplitude of the normal loading. The minimization of 


the potential energy with respect to the displacement coefficients U 


X 


and W then yields the following set of linear equations 
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For a symmetric three-layer plate with 


(4.4a) 


(4.4b) 


(4.4c) 


(4.4d) 
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A - 
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the solution of Eqs. (4.4) may be written as 
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2L 
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(4.5c) 




(4.5d) 


Direct stresses may be obtained from Eqs. (2.3) and (3.5) as 




j _ \ du 


du 


i+1 


dx 


(4.6) 


Transverse shear stresses may be obtained from the solution of the 


equilibrium equation 
3a^(C) 3Ti.(C) 


JEZ. 


3x 


35 


(4.7a) 


yielding 
i 


T (5) 
xy 


( 0 ) 

xy 


•ii= 


(- i« 


a'^ i 
d u 

2E + _S_ 

^ 2 2h. 

dx i 


,2 i+1 
d u 

X 


dx 


(4.7b) 


together with the boundary conditions of zero shear stress at one of the 
surfaces. The remaining zero shear stress conditon is then automati- 
cally satisfied. For the case of three symmetric layers 
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■ 25 X 10^ psi 
- 1 X 10^ psi 
“ 0.5 X 10^ psi 
» 0.2 X 10^ psi 

’’l " **2 * T 
i-4,10 


A comparison of the results of Ref. 13 and those given by the present 
analysis are shown in Figs. 4-2 to 4-5. It will be seen from these figures 
that despite the grossness of the approximations the midsurface dis- 
placements and stresses are in quite good agreement. The accuracy could 
be Improved even further by considering each layer to be divided into a 
number of sublayers, each of which is analyzed on the basis of the 


present theory. 

The finite element method using the triangular element outlined in 
section 2.7 was also used to solve the problem. A half-strip of the 
plate was divided into right triangular layer elements as shown in Fig. 
4r6 where the numbering of nodes and layer elements is shown. Two 
numbers shown in a triangle indicates that the interface belongs to the 
two adjacent layer elements that are numbered. Only half of the plate 
strip need be used since the plate center is an axis of symmetry with 
« 0. A state of plane strain was imposed by requiring that dis- 
placements u and w at corresponding points on both long edges of the 
strip be identical whicle the displacements u^ were taken equal to zero 
at all points. The displacements w at the top and bottom of each layer, 
while numbered differently in the layer element are made to be equal by 
assigning each the same number in the overall numbering system. 

The width of the strip was taken as 1 inch and the thickness of 
each layer was taken as 2 inches. The plate length was defined to be 
24 inches and 60 inches corresponding to length- total thickness ratios 
of 4 and 10, respectively. The material constants in reference 13 are 
given as 

Ej - 25 X 10^ psi 
E 2 ■ 1 X 10^ psi 

0^2 “ 0.2 X 10^ psi 
'’12 * ^23 


0.25 


..Jl 





Since 


'*12 ^2 " ^21 


the Poisson’s ratio is equal to 


V,, . 0.25 X ^ - ^- ", - O.OI 

25 X 10® 

SO that for layers I and 3 

25 X 10^ .. ,«6 

^11 “ 1^25' x' O.bl “ ^ 

^2 1 X 10^ , ,«6 

En-i “ ■; — :: r"A ~ Tc — “ l.oo x 10 psi 

22 ^”'^i'>'^21 1”0*25 X O.Cl ^ 


®12 ■ '' l 2®^22 ■ '’ 21^11 ■ 


E 33 ■ 0.5 X 10 psl 


» 0.5 X 10 psl 


G 33 ■ 0.2 X 10 psl 


while for layer 2 


1.0 X 10 psl 


25.06 X 10 psi 


0.25 X 10 " psi 


0.50 X 10 psi 


“ 0.2 X 10^ psi 
G^2 ■ 0.5 X 10^ psi 


ORIGINAL PAGE IS 
OP POOR quality; 


The actual half sine wave load distribution was used in determining the 
work expression. For a state of plane strain with w in each element 
given by 




the work expression for each triangle is then 


«i - i ’o 7 


m ^ 

L " x.-x, L 3 . 


IT ^2"*1 




L21S1 

L 2 


7 T ’' 2-^1 ^^^ 2 ' 

L “2“ “ . 


(4.9) 


where b is the width of the plate strip. 

Some results obtained by solvit g the finite element equations with 
the aid of the developed computer program are shown in Figs. 4-2 to 4-5 

where they are compared to the exact solution of the equations of the 
present theory. It will be seen that the finite element model used is 
capable of yielding accurate results provided a sufficiently large num- 
ber of elements is used. 

Also shown in the figures are the results obtained by Spilker, et 
al. (ref. 19), using a similar hybrid-stress finite element model which 
has the same displacement distribution but which utilizes an assumed 
stress variation in the derivation of the element stiffness matrix as 
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well. It will be seen that the present element model yields results 
that are about as accurate as the somewhat more complex hybrid-stress 
element. 

4.2 Two-Layer Angle-Ply Square Plate Under Uniform Load 

As a check on the accuracy of the present finite element model 
when applied to laminated plates, the problem of a uniformly loaded 
square plate with two angle-ply layers was Investigated. The fiber 
orlent^i^ions of the layers are +Q for one and -0 for the other. The 
layers are of equal thickness, h/2. The plate Is supported on all four 
edges, so that transverse displacements are not permitted to occur. 
Interlayer displacements normal to the edge are not permitted, but 
displacements parallel to the edge are unrestrained. Normal and parallel 
edge displacements of the top and bottom surfaces of the plate are free 
to occur. Material properties, dimensions and geometry are shown In 
Fig. 4-7. Whitney (ref. 34) has solved a similar problem with classical 
lamination theory in which shear deformations are neglected. Spilker, 
et al. (ref. 17), have solved the same problem using a simplified 
hybrid-stress element and have obtained good agreement with a corrected 
version of Whitney's results. 

Due to coupling of bending and extension In the cross-ply laminated 
plate, there are no symmetry conditions that can be applied. Thus the 
entire plate must be modeled. An isotropic single layer square plate 
under uniform pressure was investigated first. In this case symmetry 
conditions permit the modeling to be limited to only one eighth of the 
plate. For a plate with a length-thickness ratio of 10, the finite 
element results converge to a value about 6% greater than the value 
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given in ref. 31 (Fig. 4-8) when each half-length was divided into 
20 equal parts, i.e. 100 elements were used for 1/8 the plate. The 
larger value to which the finite element solution converges is due to 
shear deformations which are not Included in ref. 31. 

For the cross-ply laminated square plate, the above results indi- 
cate that each side of the plate should be divided into 40 equal parts 
to achieve satisfactory convergence. Thus 3200 layer elements would be 
needed to model the entire plate. Results that were obtained with the 
side divided into 6, 8 and 10 spaces respectively (Fig. 4-9) are in 
accord with this conclusion. The hybrid-stress element developed by 
Spilker, et al. (ref. 17), provided convergent results for the angle- 
ply laminated plate with the side divided into 10 spaces. It is 
apparent, then, that a plate finite element model based upon the prin- 
cipal of m-t nimiim potential energy has very much poorer convergence 
properties. 
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4.3 Perforated Square Plate Under Uniform Inplane All-Round Tension 
The effect of circular holes on the stress distribution In 
stretched plates has been the subject of many Investigations. Accurate 
solutions are avallabe for a circular hol4^1n an Infinite plate 
(refs. 34-37) and In an Infinite strip (refs. 38-39). The effect of 
circular holes In a square plate has been studied less rigorously by 
Hengst (ref. 40) vdio assumed an expression for the Invariant trace of 
the stress tensor of the form 

4 

P^r^^ cos 2k9 
k-4' 

and obtained expressions for the coefficients by minimizing the error of 
the expression. The investigation yields an approximate formula for 
the circumferential stress around the hole. For uniform tension in 
both directions (Fig. 4-10) the circumferential stress for various 
size holes is given in Table 4-1 . 

For the finite element solution of a plate with r^/a » 0.5 the 
plate was modeled as shown in Fig. 4-11. The symmetry of the problem 
requires that only one-eighth of the plate be modeled. The displace- 
ment symmetry conditions for the problem are those of identical dis- 
placements in the x- and y-directions along the diagonal a-b and van- 
ishing y-direction displacements along the line c-d. Along the side 
b-c and the circle a-d the stresses are prescribed. The symmetry con- 
dition of identical displacements is achieved by assigning the same 

variable number to U and U at a nodal point on the diagonal. The 

X y 

problem is one of generalized plane-stress with stress boundary condi- 
tions so that the results are independent of material properties and of 
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accuracy of the model since the prescribed and calculated stress values 
need not be Identical. The circumferential stress values along the 
centerline are Interesting In that they decrease at first but then 
Increase somewhat toward the loaded edge. 


4.4 Perforated Square Plate Under Uniform Nodal Pressure 

The square isotropic plate with a circular hole which Is subjected 

•1 

to uniform load has been Investigated by El-Hashlmy In reference 41 
A deflection function Is assumed of the form 


w 




+ (A^p ^ + B^p ^ + B 2 P^) cos 40 


+ (A 3 p“® + B 3 p"^ + A^p® + cos 80] 

where p Is the non-dimensional radius ratio r/a. The hole radius Is 

p^a. Here a is the side length of the plate. The 12 coefficients are 

determined by the satisfaction of the conditions of vanishing moment 

and effective shear force at the hole, which yields 6 equations, and 

by the satisfaction of conditions of simply supported or clamped edges 

at 3 points as indicated in Fig. 4-15. The periodicity of the solution 

requires that only an eighth of the plate be considered. 

The finite element method of the present paper was used to analyze 

the plate as well. The hole size and the grid used are identical with 

those of the previous section. The boundary conditions are those of 

zero force at points on the hole, equal values of U and U and zero 

X y 

effective transverse shear force along the diagonal, zero values of U^ 
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and effective transverse shear force along the centerline, and zero 
transverse displacement and either zero Inplane stresses displacements 
along the edge, corresponding to simply supported and clamped edges 
respectively. The number of degrees of freedom are then 648 for the 
simply supported plate and 376 for the clamped plate. 

A comparison of the transverse displacements calculated using E1-* 
Hashlmy's equations and those obtained from the finite element anlaysis 
Is shown in Table 4-2. It will be seen that the finite element values 
are considerably less than those of ref. although the shapes of the 
deflection patterns are similar. This result Is attributed to the sat- 
isfaction of boundary conditions at only 3 points on the edge in El- 
Hashlmy's solution which yields a much more flexible structure than is 
actually the case. 

The deflections and stresses obtained from the finite element 
method are shown in Figures 4-16 to 4-17 • 

4.3 AxLsymmetric Bending of a Cylindrical Shell by End Load 

Having ascertained by the investigation of the previous problems 
that the results for flat plates obtained by means of the present 
finite element model are adequate, the accuracy of the theory for 
cylindrical shells was studied next. The first problem investigated 
was that of axisjnnmetric bending of an isotropic cylindrical shell by 
end load, with transverse shear deformations included. The general 
equations of the theoEy for a shell having ortho tropic layers with 
the material axes coincident with the lines of curvature axes are given 
in Appendix C together with the solution of these equations for a 
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single layer shell loaded equally at both ends by (?lg. 4-18^) radial 
shear forces. 

Since the deformations are Independent of position around the 
cylinder, elements of any width may be used In a finite element solu- 
tion. The central symmetry of the problem allows the modeling to be 
limited to half of the shell to be modeled. Thus a portion of the 
shells of length equal to one-half the total length and a width of 1 
Inch was divided Into finite elements as shown In Fig. 4-19. To study 
the convergence of the procedure the nxmiber and size of the elements 
was varied. In some cases the size of the elements was equal. In 
other cases a portion of the shell nearer the load was divided into 
p smaller elements while the remaining portion near the center was 
divided Into m larger elements. The maximum number of equations to be 
solved was limited to about 500 in the present calcxilatlons . The mate- 
rial of the shell was taken to be transversely isotropic with Young’s 
modulus E equal to 1000 ksi and Poisson's ratio equal to 0.3. The 
transverse shear modulus was taken as 100 ksi. The thickness of the 
shell was assumed to be 1 inches. The inner radius was varied from 

5 inches to 100 inches. The shell length varied with the shell thick- 

L_ 

ness so as to maintain the value of about 6.5. in 

the calculations axisymmetry was imposed by reqxxiring that the dis- 
placement Urt be equal to zero everywhere and that u and w be indepen- 
tj x 

dent of circumferential position. The boundary condition of symmetry 
about the center of the cylinder was Imposed by making, the displacement 
and the transverse shear force equal to zero. At the loaded edge 
the axial stress was put eq\ial to zero while the radial shear force 
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was kept constant for all shells. The results of the calculations for 
the end dtvflectlon are given In Table 4-3 together with the values 
calculated for the exact theory of Appendix C. 

The rate of convergence of the finite element results is shown in 

Fig, 4-19. where the ratio of the finite element and anlaytic maximum 

# 

displacements is plotted as a function of the number of elements in the 
half-length of the shell. It should be noted that the number of ele- 
ments shown is that of the triangular elements which is thiis twice the 
number of linear a3d.symmetric cylinder elements. The results indicate 
that the thicker shells require less elements to achieve accurate dis- 
placements but that sufficient accuracy can be obtained within the 
equation limitations set. 

The displacements and direct stresses obtained with the finite 
element analysis xising the largest number elements are compared with 
analytic results for a long cylinder in Figures 4-21 to 4-26. The 
stresses were obtained by calculating the derivative of longitudinal 
displacements as for the beam element and combining this with the aver- 
age radial displacement within the element. The very good agreement 
between finite element and analytic results indicates that a relatively 
unsophisticated element model will yield satisfactory results. 

4.6 Cylindrically Curved Square Plate With a Circular Hole Under 

Uniform Transverse Pressure 

A final problem to be treated is the effect of curvature on the 
deflections and stresses of a single layer square plate with a central 
circular hole. The developed geometry and material properties of the 
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curved plate are similar to those of the flat plate of Section 4.4 
except that the plate has a cylindrical radius of 1000 Inches. Since 
the geometry parameter z* (*^i2(l-\)^) ^ s m all, the results should 

be similar to those for the flat plate. 

One quarter of the curved plate was modeled as shown in Fig. 4-27 
since the diagonal is no longer an axis of synmetry. A comparison of 
the radial deflections for the curved plate and those for the flat 
plate is given in Table 4- 4 . The deflections are comparable but even 
a slight amount of curvature is sufficient to change the deflected 
shape significantly and to reduce the maximum deflection by 13%. 

The distribution of stresses at the inner and outer surfaces of 
the curved plate are shown in Figs. 4-28 to 4-29. 



CHAPTER V 


CONCLUSIONS AND RECOMMENDATIONS FOR 
FUTURE INVESTIGATION 


The results of the present investigation have indicated that the 
addition of shear deformations to the theory of laminated plates and 
shells permits the use of a simple linear displacement variation finite 
element model for calculations. Provided enough elements are used, 
accurate displacements and direct and shear stresses parallel to the 
layer surfaces can be obtained. Where conditions of symmetry permit 
the analysis to be limited to that of a small region of the structure, 
the elements can be made small enough while the total number of equa- 
tions is kept reasonably small. However, for larger structures, the 
use of such a simple model in conjunction with the principle of minimum 
potential energy sacrifices rapidity of convergence for simplicity of 
calculations as compared to a similar hybrid— stress finite element 
model due to Mau, Tong and Plan. A linear displacement minimum poten- 
tial energy element yields strains and stresses which are constant in 
each element and hence it does not permit the direct calculation of 
transverse shear and normal stresses which depend on the derivatives of 
the direct and shear stresses parallel to the layer surfaces. The 
hybrid-stress element, on the other hand, does permit the direct calcu- 
lation of trnasverse shear stresses. 
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One would be tempted to abandon the potential energy approach were 
it not for the ease with wlch the approach Is extended to shell analy- 
sis. The extension of the hybrid-stress model Is difficult, requiring 
the definition of a stress state within the shell wall which satisfies 
the equation of equilibrium and Interlayer continuity conditions. 

Such equilibrium stress state distributions would vary with the shell 
geometry and would be difficult to define. Other disadvantages of the 
hybrid-stress model are the large amount of computer storage required 
to define the stress-state for each element as well as the numerous 
matrix Inversion calculations required to obtain the combined stiffness 
matrix for a stack of layer elements extending through the wall thick- 
ness. These calculations Increase with the number of layers. A topic 
for future Investigation that should be pursued, therefore. Is the 
effect of more complicated displacement variations on the accuracy and 
rapidity of convergence of potential energy finite element models. 

The right triangle element chosen in the present investigation for 
ease of derivation and calculation is useful for regular structures. 

It Is inconvenient for irregular structures, however, since it does 
not permit complete flexibility in the use of small triangles in a 
region of rapid stress variation and larger triangles in regions of slow 
variation. With the right triangle, the small triangle dimensions 
affect the size of the larger triangles. Thus the analysis should be 
extended to triangles of arbitrary shape. 

In the absence of these extensions, the analysis technique de- 
veloped in the present paper is useful but may require the solution of 
an extremely large number of simultaneous equations for accurate 
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dlsplacemeivt: and stress states in layered structures. 
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Table 4-1 Variation with Hole Size of the Circumferential 
[1 Stress for a Perforated Square Plate Under Uniform 

! Tension 




09/q 


1/8 


2.028498 - 0.000707 cos 40 - 0.000000 cos 86 


1/4 


2.119101 - 0.011014 cos 46 - 0.000010 cos 86 


3/8 


2.289766 - 0.053591 cos 40 - 0.000237 cos 86 


1/2 


2.582663 - 0.164303 cos 40 - 0.002353 cos 86 






/i 


Finite 

Ref. 41 Element 
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10.692 


3.679 


0.000 



85.498 


69.414 


33.876 


0.000 


5.569 


3.293 


1.247 


0.000 


Finite 
Ref. 41 Element 



15.171 


11.889 


5.769 


1.264 


0.000 


87.242 


72.581 


44.113 


17.518 


0.046 


5.735 


4.015 


1.867 
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39.483 


33.265 


17.952 


6.300 


Finite 

Ref. 41 Element 


15.438 


12.299 


6.567 


2.256 


0.450 

0.168 


87.825 


73.644 


46.650 


23.918 


9.934 


0.000 


6.227 


5.200 


2.635 



39.873 


33.990 


21.831 


11.156 


6.877 


0.000 


Table 4-2 Comparison of the Transverse Displacements of a Perforated Square Plate Under 
Uniform Transverse Pressure (rn/a ■■ 0.5) 
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Table 4-2 Comparison of the Transverse Displacements of a Perforated Square Plate Under 
Uniform Transverse Pressure (x^/a = 0.5) 


























Nuabtr of| t-U.O", | L-16.0", »-IO. 

Olvlalona 


L-2J.O", 1-20.0" I L-2fc.O", »-50.0" 
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Table 4-3 Comparison of Finite Element and Analytic Values of End Deflection of a (^lindrlcal 
Shell under Axisymmetrlc End Load 




























Comparison of Curved and Flat Plate Normal Deflections 
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C.P. - curved plate 
F.P. - flat plate 


Table 4-4 Comparison of Kadlal Deflections for Perforated Flat and Curved Plates Under Uniform 
Pressure 
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Figure 2-7 

Numbering system for triangle vertices. 



Figure 2-8 

N\moberxng system for vertex displacements. 
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Figure 2-9 Different Element Orientations 
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GEOMETRIC PROPERTIES: a-b-10" 

MATERIAL PROPERTIES; E^^-40xl0®psi E^2-l x lO^si 

Gj^ 2 "<' 23 "®*^ * lO^si 

V12-V23-O.25 

FIBER ORIENTATIONS: -6 In layer 1; +6 in layer 2 

LOADING: Uniform transverse load, 100 psl 

BOUNDARY CONDITIONS: Simply supported on all four sides 

FINITE ELEMENT MESH: MxM square elements in entire plate 


Figure 4-7 Laminated Plate Under Uniform Transverse Pressure 
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Figure 4-16 Deflection Distribution for a Perforated Square 
Plate with Clamped Edges 
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Figure 4>20a Convergence R^te of Finite Element Solution 
for End Displacement 
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Figure 4-25b Comparison of Finite Element and Analytic 

Stress Distribution for a Shell with R/h"50 
Bottom Surface 



Figure 4-26a Comparison of Finite Element and Analytic 

Stress Distribution for a Shell with R/h»100 
Top Surface 
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Figure 4-26b 


Comparison of Finite Element and Analytic 
Stress Distribution for a Shell with R/h*I00 
Bottom Surface 
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Figure 4-27 Finite Element Model for Curved Plate With a 
ClrcTilar Hole 
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iMjfure 4>28b Inner Stress a« Distribution for Curved Plate 

°o 

with a Circular Hole 
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Figure 4-28c Inner Stress Distribution for Cur:.ed Plate 

o 

with a Circular Hole 






















































Figure 4-29a Outer Stress a Distribution for Curved Plate 
with a Circular Hole 
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Figure 4-29c Outer Stress Distribution* for Curved Plate 

with a Circular Hole 
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